R包poLCA的潜类别分析（LCA）
poLCA的LCA
先看一个简单示例。
library(poLCA)#数据集，详情 ?values
#调查了 216 名受访者对四个问题（A、B、C、D）的认同倾向，“1”代表不认同，“2”代表认同
#注：poLCA 中，分类变量使用 1、2、3 等数值代替，不同数值表示不同类别
data(values)
head(values)
#执行 LCA 识别数据集变量间的潜类别结构，详情 ?poLCA
#例如设置潜类别为 2
#从随机起始值开始，1000 次迭代获取参数估计值（对数似然值）
#并运行 5 次独立的 poLCA 获得可能的全局最大值
M1 <- poLCA(formula = cbind(A, B, C, D)~1, data = values, nclass = 2, maxiter = 1000, nrep = 5, graph = TRUE)
M1
poLCA()的formula书写为cbind(dv1,dv2,dv3...)的样式，dv1、dv2等代表观测分类变量。对于没有协变量的情况，公式为cbind(dv1,dv2,dv3)~1；若考虑协变量，将~1替换为所需的预测变量函数，例如cbind(dv1,dv2,dv3)~iv1。
poLCA()使用局部独立性假设估计潜在多向表（latent multi-way tables）的混合模式，其数量（nclass）需手动指定。估计的参数包括每个外显变量的条件概率，“混合”比例表示与每个潜在多向表相对应的观测值的population share，以及在公式中指定的任何分类预测协变量的系数（存在协变量时）。
poLCA()使用期望最大化（Expectation Maximization，EM）和Newton-Raphson算法获得潜类别的最大对数似然（log-likelihood）参数估计值。受起始值（参数probs.start默认NULL，每次计算将使用不同的随机初始值）和迭代次数（参数maxiter默认1000，即迭代1000次）的影响，它可能仅定位局部最大值，而非全局最大值，并随nclass设定值的增加更可能出现问题。因此建议多次运行poLCA，直到相对确定找到了全局最大值的可能性，例如将参数nrep的值设置为大于1（默认1），可以通过一次计算多次估计LCA的参数值并自动选择全局最大值。
条件概率和类别概率
结果概要如下所示。
首先通过5次计算结果的参数估计值可知，已经获得了全局最大值，所获得的两个潜类别是有代表性的。
Conditional item response probabilities中，列出了各分类变量中各外显变量在潜类别中的条件概率（每一行的概率之和为1），可据此判断哪些潜类别主要与哪些外显变量有关。
Estimated class population shares，为本示例所获得的两个潜类别的估计概率，可以将潜类别概率可近似理解为因子分析中的解释方差比例，以评估潜类别的重要程度。
通过拟合所获得的潜类别至原数据，获得AIC、BIC等适配检验指标，可以用于帮助评估潜类别是否合适。例如，分别比较nclass=n和nclass=m（n<m）的结果，综合考虑AIC、BIC值更低且潜类别更少的情形。
summary(M1) #对象中的各元素所记录的内容，详情 ?poLCA
M1$probs #Conditional item response probabilities
M1$P #Estimated class population shares
#即 Conditional item response probabilities 中的值
library(reshape2)
library(ggplot2)
M1_probs <- melt(M1$probs, level = 2)
ggplot(M1_probs,aes(x = L2, y = value, fill = Var2)) +
geom_bar(stat = 'identity', position = 'stack', width = 0.5) +
facet_grid(Var1~.) +
scale_fill_brewer(type = 'seq', palette = 'Greys') +
theme_bw() +
labs(x = '', fill = 'probabilities') +
guides(fill = guide_legend(reverse = TRUE))
后验类别概率
潜类别分析还可用于实现对象聚类。
在找到理想的潜类别并获得参数估计后，继续使用上述两类概率值（外显变量的条件概率和潜类别概率）计算外显变量反应每种组合分类到各潜类别的后验概率，然后根据后验概率的大小决定该组合应归入的潜类别，即创造一个新类别变量说明观察值的后验类别属性。
#后验类别概率M1$posterior
#对象归类
M1$predclass
#作图
library(ggplot2)
library(ggrepel)
posterior <- data.frame(M1$posterior)
posterior$label <- rownames(values)
posterior$class <- as.character(M1$predclass)
names(posterior)[1:2] <- c('class1_probabilities', 'class2_probabilities')
ggplot(posterior, aes(class1_probabilities, class2_probabilities, color = class)) +
geom_point() +
geom_text_repel(aes(label = label, color = class), size = 3, show.legend = FALSE)
结果用于表示对象属于某一类别的概率，并根据对象后验概率的最大值将其归入特定类别。例如上述显示对象1（对应于原数据集，即第一位受试者）在class2中的后验概率最大，因此将它归类为class2类，其它同理。
带协变量的LCA
以下再展示一例考虑协变量时的LCA示例。
带协变量时，LCA也可称为“带协变量的潜类别分析”（latent class analysis with covariates），或“伴随变量潜类别分析”（concomitant-variable latent class analysis），其中潜类别的可能性由一个或多个协变量预测。
## Three-class model with a single covariate.data(election)
nes2a <- poLCA(formula = cbind(MORALG, CARESG, KNOWG, LEADG, DISHONG, INTELG, MORALB,
CARESB, KNOWB, LEADB, DISHONB, INTELB) ~ PARTY, data = election, nclass = 3, nrep = 5)
nes2a
pidmat <- cbind(1, c(1:7))
exb <- exp(pidmat %*% nes2a$coeff)
matplot(c(1:7),(cbind(1,exb)/(1+rowSums(exb))), ylim = c(0,1), type = 'l',
main = 'Party ID as a predictor of candidate affinity class',
xlab = 'Party ID: strong Democratic (1) to strong Republican (7)',
ylab = 'Probability of latent class membership',
lwd = 2, col = 1)
text(5.9, 0.35, 'Other')
text(5.4, 0.7, 'Bush affinity')
text(1.8, 0.6, 'Gore affinity')
此时在poLCA()的formula中，将cbind(dv1,dv2,dv3)~1（无协变量情况）替换为cbind(dv1,dv2,dv3)~iv1（协变量iv1）之类的形式。更复杂的其它情况，如cbind(dv1,dv2,dv3)~iv1+iv2*iv3，考虑iv2和iv3交互作用。
与上述相比，结果输出中还给出了类别预测协变量的系数等值，可帮助进一步评估协变量对预测潜在类别成员的影响。